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We develop analytic solutions for the linear evolution of metric perturbations in the DGP 
braneworld modified gravity scenario including near-horizon and superhorizon modes where so- 
lutions in the bulk are required. These solutions apply to both the self-accelerating and normal 
branch and elucidate the nature of coordinate singularities and initial data in the bulk as well as 
their effect on perturbation evolution on the brane. Even on superhorizon scales, the evolution of 
metric perturbations is no longer necessarily scale free due to multiple resonances in the bulk. Based 
on these analytic solutions, we devise convenient fitting functions for the evolution that bridge the 
^. . various spatial and temporal regimes. Compared with a direct numerical integration of the bulk 

equations, the fits are accurate at the percent level and are sufficient for current and upcoming 
observational tests. 
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I. INTRODUCTION 



£""N | The Dvali-Gabadadze-Porrati (DGP) model modifies General Relativity on large scales by positing that we live on 

r \ ■ a 4-dimensional brane in a 5-dimcnsional Minkowski bulk |l|. O n scales larger than the crossover scale r c , gravity 

becomes 5-dimcnsional and hence one can hope to uncover extra-dimensional physics by studying the evolution of 

^G . density perturbations on scales approaching r c . The DGP model has two branches of cosmological solutions [3|. On 

Mh' the self-accelerating branch, the expansion of the universe accelerates without a cosmological constant or dark energy. 



q , On the normal branch, brane tension is required to accelerate the expansion though gravity is still modified on large 
?— i ■ scales. 

To solve for the evolution of perturbations on either branch near the horizon scale and beyond requires following 

Cw ■ metric perturbations into the bulk. On the self-accelerating branch approximate scaling solutions Q showed that 

modifications near the horizon scale are in significant tension with the data from the cosmic microwave background 

CNJ ' tallf- I n this paper, we present the first quantitative and thorough comparison of these approximate solutions with 

£> ' direct numeric simulations [6[ . We confirm that the scaling solutions used in [J, Q are sufficiently accurate to expose 

O^l tension between self-accelerating DGP and observations. It is also worthwhile noting that the pure de Sitter phase 

of the self-accelerating branch suffers from a perturbative ghost mode (see e.g. Q), which is significant theoretical 

^ I | challenge for the model. 

. ■ Numerical integration of the bulk metric equations have also been performed on the ghost-free normal branch 

Q. In order to facilitate studies of cosmological constraints on the normal branch, it is useful to have a simple 



but accurate description of the evolution of perturbations in terms of the fundamental cosmological parameters and 
crossover scale. Scaling approaches on the normal branch have also been studied [8[ but without a proper treatment 
of initial conditions in the bulk and boundary conditions on the brane. 

In this paper, we employ analytic and numerical techniques to devise such a description. In fact, preliminary 
versions of results from this work have already been used to set cosmological constraints on the crossover scale on 
both branches of DGP @. 

We begin in JT] with a brief review of the equations governing the background and perturbations on both the self- 
accelerating and normal branches of the DGP model. We outline the numerical methodology for obtaining solutions 
of the perturbation equations as well as the scaling ansatz that provides insight into their nature in ijllll In ^IVI 
we combine scaling results with Green's function techniques to obtain analytic solutions in the matter and de Sitter 
epochs. We join these solutions into simple but accurate global descriptions of perturbation evolution via fitting 
functions in ffV] and discuss these results in WII 

II. FORMALISM 

We begin in ^11 Al with a review of the DGP field equations and bulk geometry. In mi Btill Cl wc apply the field 
equations to the background evolution and perturbation evolution respectively. Wc discuss the initial conditions and 
singularities in the bulk in ^11 Dl and the effective equations of motion on the brane in ^11 El 



A. Field Equations and Bulk Geometry 



The action of the DGP model is 



S = 7^2 J <PXj—gR^) + _L j d *x^R {i) + J d A x^j(C m - a). (1) 

M dM b dM b 

Here, a is the brane tension and L m is the matter Lagrangian. The field equations satisfied in the empty bulk are 
simply 

R'S = 0, (2) 

whereas those on the brane are given by [1JJ, [HI 

G$ = (2^ rc ) 2 n^ - S„ v , (3) 

where 

T^v =T >1V - ag a p - k± Gffi, (4) 

and £ M „ is the trace-free projection of the 5-dimensional Wcyl tensor. Here the crossover scale is defined by 

- - k < 5 » 

The background metric or geometry converts the field equations into the modified Fricdmann equation and imposes 
the distinction between the two branches of solutions. The background bulk metric can be expressed in terms of 
Gaussian normal (GN) coordinates (t, y, x) [12[ 

ds 2 = -?i 2 (i, y) dt 2 + 6 2 (f, y) dx 2 + dy 2 , (6) 

or null coordinates (u, u,x), 

ds 2 = -du dv + r~ 2 v 2 dx 2 . (7) 

For the analytic solutions below, we employ GN coordinates whereas for the numerical characteristic integration (CI) 
we use null coordinates. 

In GN coordinates the brane is always at y — 0, while in the null coordinates the brane is moving. The metric 
functions in the GN line element are 

n(t,y) = l + e(H + ^)y, b(t,y) = a(t)(l + eHy) , (8) 

where e = +1 on the self- accelerating branch and e = — 1 on the normal branch. Here H = a /a is the Hubble 
parameter on the brane. 

The coordinate transformation between the two systems is 



1 

u = — 



a(t) 



da ey 



H 2 (a)a? Ha 



+ U , v = r c b(t,y), (9) 



where u(ao) = uq is a constant defining the origin of the coordinate system at ao- Ignoring inflation and taking only 
matter components on the brane, we can take ao — > and set uq = 0. 

With these assumptions the coordinate transformation fixes the position of the brane in the null system: 



Ub(t) = — I 
r c Jo 



1 f a(t) da 

H 2 {a)~a 2 



v h (t)=r c a(t). (10) 



Note that both line elements have a coordinate singularity: either when b(t, y) — or v — 0. These singularities are 
actually at the same place. The relationship between the two coordinate systems is illustrated in Fig. Q] Note that 
the v = singularity is only accessible in the normal branch. 
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FIG. 1: Spacetime diagram of a brane trajectory through the bulk for the self-accelerating and normal branch. The grey 
timelike lines are the y = constant traces of the Gaussian normal coordinate system, while the spacelike grey lines are the 
t = constant traces. Lines of constant u and v are tilted at 45° and represent the trajectory of light rays in the bulk. The 
brane is at y = 0. In both cases there is a past Cauchy horizon H- in the bulk, and for the normal branch %_ is coincident 
with a coordinate singularity. In this plot, we have chosen the constant «o in Eq. ((9| such that the Cauchy horizon in the 
self-accelerating branch is at u = 0. 



B. Background Evolution 

The field equations in the spatially flat background reduce to 

1 



H = 



2r c 



e + \/ 1 + 3^r c 2 (p + <J) 



(11) 



with the bulk metric determined by H through Eqs. (|5HTu]). If we further assume that the energy content of brane p 
is dominated by non-relativistic matter, this modified Fricdmann equation can be rewritten as 



tr 

— - = \/f2 m a- 3 + Oa + VL Tc + e\/Tlj~, 

where 51a = kI<t/3Hq is the effective cosmological constant for the brane tension and 

1 e 



'Sir 



-(1 — VL m — 0/ 



2H r c 2 
defines the crossover scale in terms of the other cosmological parameters. 



(12) 



(13) 



C. Perturbation Equations 



The perturbed GN line element is written as: 



ds 2 = -n 2 (l + 2XA) dt 2 + 2n\A y dtdy + b 2 (l + 2XR) dx 2 + (1 + 2\A yy )dy 2 



(14) 



Ay = 

A vv = +^T ( "HTT + ~~^PT ~ ~~~H7 ~ r~~ ) ' ( 15c ) 



4 

in the longitudinal gauge where A <C 1 is an order parameter for bookkeeping purposes. These metric perturbations 
are governed by the master variable Q [13[ through the relations [lj] 

1 ( 9 d 2 n 1 d 2 Q 1 dndn ldndft\ 

66 \ dy 2 n 2 dt 2 n 3 dt dt n dy dy J 

1 / d 2 n ldndQ. 
bn \ dtdy n dy dt 

1 fd 2 n 2 d 2 fl 2 dndfl 2dndfl 

66 \ dy 2 n 2 dt 2 n 3 dt dt n dy dy 

1 fd 2 n 1 d 2 fl 1 dndfl ldndfl\ 

66 \ dy 2 n 2 dt 2 n 3 dt dt n dy dy J 

The master variable satisfies the wave, or master, equation 

d ( 1 d®\ d fndn\ n l2n „ ,„ . 

-m{^m) + d- y {vl) y -)-¥ kn = > (16a) 

d 2 n A^_^l^-o (16b) 

dudv 2v dv 4v 2 

in GN and null coordinates, respectively. 

The master variable also satisfies the boundary condition 



£2lo 4. ^Hlo - 3 ( e ^fc 2 +74g 2 a 2 ) 3er c ^pa 3 74 

2H b 4 b AHa 2 b 2k 2 



(C y S2)b--7 TI 7S2b + — ;— S2b 7773 S«b + ^ A - I 17 ) 



The boundary value of the master variable £\ acts as a source to the brane metric perturbations \& = A(y = 0) and 
$=ft(y = 0): 

* - + ^fc^ A + 4o^ b 12Hr c a 3 ^ (18&) 

* = -^g^A + -^ b - ^O b + ^(^c74 + ^ 72 ) 

2fc 2 4ffr r a 4a 4r r a 3 



where the dimcnsionless 7-factors 



2eHr c 

^ = T77 7' 19a 

2eHr c — 1 



2er c (H -H 2 + 2eH 3 r c ) 
H(2eHr c - l) 2 



7 2 = %;,„,;;_ ::; c \ ^ 



4er c (2er c H - 3H + 6eH 2 r c ) . 

73 = 9(2eHr c - I) 2 ' (19c) 

Ae{er c H-H + 2eH 2 r c ) 

74 = 3H(2eHr c I) 2 ' (1M) 

With these relations, energy-momentum conservation of the matter on the brane simplify to yield 



1 2 e-y 4 k 4 



A + 2i/A-- K 2 /072 A = -4V^- ( 2 °) 



Eq. ([i~6|) . (p~7)) and (|20[) represent the complete equations of motion for the perturbations. 

D. Initial Conditions 

In order to fully specify the dynamics of perturbations, the equations of motion must be augmented by initial 
conditions, not only on the brane but also in the bulk. On the brane, the suppression of DGP modifications at 
Hr c 3> 1 means that initial conditions can be set as usual with adiabatic conditions. In the bulk, it is physically 



sensible that all bulk perturbations are causally generated by perturbations on the brane; i.e., there should be no 
sources for $7 except the brane itself. This is equivalent to demanding that the bulk master variable vanish on the 
past Cauchy horizon: 



fl(«_)= 0. 



(21) 



While the principal motivation for this assumption is causality, it is interesting to note that in the normal branch this 
is a necessary condition to avoid a curvature singularity in the bulk. To see this we substitute Eqs. (|15[) and (|16a[) 
into (|14| . and then calculate the full Ricmann tensor in the normal branch. We find that the Ricmann tensor tends 
to diverge on the past Cauchy horizon H_ (where fo = and v = 0) in the normal branch unless f2 is correspondingly 
suppressed; for example, 



R 



tyty 



3 b 5 



0{\ 2 ). 



(22) 



We shall see that this issue is related to an irregular singular point at the horizon at finite k in the master equation 
and requires that Q, be exponentially suppressed near H-, which is actually a stronger condition than Eq. (I2T1) . 



E. Brane Gradient and Metric 



The dynamical equations (fTT|) and (|2T)]) are closed on the brane if the relationship between (d y il)b and fib is known. 
From the perspective of the observer on the brane, the wave equation (|16p for the master variable in the bulk simply 
serves to specify the dimensionless gradient 



R 



\eHnJ, 



(23) 



We shall see that scaling and Green's function arguments can be used to extract this gradient in various limits. 

Ultimately, we are most interested not in the brane gradient i?, but the metric perturbations themselves. In 
particular the ratio of metric perturbations 



g(a,k) = 



$ + * 



(24) 



plays a central role in modified gravity models as it vanishes in General Relativity if there is no anisotropic stress. 
Moreover in the parameterized post-Friedmann (PPF) formalism |15| . it efficiently specifies metric perturbations for all 
scales. For the large scale regime (k/Ha) <C 1, in General Relativity it is well known that the curvature perturbation 
on comoving slices £ is conserved. This conservation law is a consequence of energy-momentum conservation and 
applies to DGP as well Q. For the DGP model, we can express ( in terms of fl^ and A: 



c = $ 

The derivative of £ can then be shown to be 



Ha, 

6q 

P 



Ha 2 



-A 



Ha 2 

k 2 



-A 



£7i fc2 
12Hr c a 3 



a, 



S 12 V Ha 



Hr c 



7i f2 b + (71 - 97 3 )O b 



(25) 



(26) 



where a prime indicates differentiation with respect to In a. One can also confirm that the following identity holds in 
DGP: 



H 1 \ H W ) 12 V Ha 



Hr c 

a 



ng 



f2 b + 2H'H- 2 { 11 H I + 2H)Qi 



(27) 



Hence, we recover that 



in the k ->• limit @. Given g, $ = ^(g 
and H. 



Till / TTl TTll \ 

C' w «$"-*' - — *' - — - — U, (28) 

H' \H H> J ' v ' 

\)/{g — 1) and this ODE can be readily integrated to find $ or \P given g 



Conversely at subhorizon scales k/Ha ~^> 1, one can employ a "quasistatic ansatz" which assumes that gradients 
with respect to x and y are of the same order and dominate over time derivatives (see ^IV Al for more details on this 
regime). This implies that the influence of the bulk through the brane gradient R is negligible, as it involves only one 
spatial derivative, and [16[ 

k*n h w 2Hr c — K 2 4 pa 5 A. (29) 

73 

Eq. (fT5)l can be expressed as the Poisson equation [17| 

V * &*■ « 

for the quantity ($ — v E')/2 which enters into observables such as gravitational redshift and lensing. The quantity 'J, 
which enters into the motion of non-relativistic matter is then specified by the relation 

a « oqs = : . (31) 

1 3[l-2Hr c e(l + H/3H 2 )} 

In both limits the dynamics of the perturbations are completely specified once g(a, k) is known. Since g(a, k) is 
determined in the quasistatic limit by the structure of the equations themselves, we concentrate on understanding the 
superhorizon regime in £11111 and mV\ We then seek a suitable interpolation between the two regimes in $Vj This can 
then be used in the PPF formalism with publicly available codes |18| to generate the metric perturbations $ and Sf! 
in an efficient manner. 

III. METHODOLOGY 

In this section we review two techniques, characteristic integration and the scaling ansatz, that have been used 
in the literature to solve numerically for DGP perturbations on large scales where bulk effects are important. In 
the following sections, we will use the scaling ansatz to develop the analytic approximations and the characteristic 
integration algorithm to test them. 

A. Characteristic Integration 

In order to model the behavior of perturbations in the DGP model, the characteristic integration (CI) algorithm 
[q.[l9l| constructs a direct finite difference solution to the bulk (fPoj) and brane (|20|) equations of motion subject to the 
boundary condition (|17[) . The code makes use of the null coordinates (u, v), which implies that the brane is a moving 
boundary. 

The natural computational domains of the CI algorithm for both branches are indicated in Figure [1] Initial data 
for the bulk field fl must be supplied on the past null boundary of the domain. In addition, one must also specify 
the value of A, A and (lb at the intersection of the initial data surface and the brane. Because of the fact that this 
algorithm is based on null curves in the bulk, the initial condition (|2ip is very easy to approximate, we merely need 
to set 51 = on the initial null surface in the computational domain. As the initial scale factor of the simulation aj n j t 
is pushed further into the past, the output of the CI code will approach the desired solution with Q('H_) = 0, though 
in practice one finds that simulation results are stable to changes in ai n i t provided that it is below some threshold 
value (typically 10 -3 for the simulations presented in the paper). 

B. Scaling Ansatz 

During epochs when the expansion rate on the brane is dominated by a single component, we expect the master 
variable on the brane to reach a scaling solution J7b oc a p for modes that are outside the horizon (k/Ha) <C 1. The 
boundary equation (fTTj) then implies that the gradient term of Eq. (|23[) obeys 

3e 3 K4P ,, A 

473P- T 74+2^ C — J.a—, (32) 




which gives the relationship between J7b and A and closes the dynamics on the branc. 

To determine the gradient parameter R we make a similar scaling ansatz for the dependence of $7 on the GN 
coordinates in the bulk Q 

£l(t,y)=a p (t)G{eH{t)y). (33) 

This ansatz is motivated by the distance in the bulk connected by a null worldlinc dy = ndt, i.e. a "horizon" in the 
bulk. For the self-accelerating branch 

y hol = aHj -^, (e = +l), (34) 

and Eq. ([9]) shows that y = j/h or corresponds to the null line u = 0. With a power-law behavior for H 2 ex a~ 3 ( 1+w ^ 2 , 
Hy hor = 1/(2 + 3w) if w > -2/3. 

For the normal branch, independently of the evolution of H 

2/hor = -^, (e = -l), (35) 

which corrects Eq. (21) in Rcf. [8|]. Eq. ([9]) shows that y = yi loI corresponds to the null line v = on the normal 
branch and hence coincides with the coordinate singularity (see discussion below Eq. (fT0|)V 

The final condition introduced in Ref. [3j is that the initial data G{eHy\ lm ) = based on the assumption that bulk 
perturbations are generated causally from brane perturbations. This is entirely equivalent to Eq. (|2ip and compatible 
with the CI initial data as discussed above. Recall that G(eiJ j/hor) = is also required to keep the Riemann curvature 
finite in the normal branch. 

The scaling technique can be extended to apply between different scaling regimes by iteratively solving for a time 
dependent p(a) with the so-called dynamical scaling (DS) method |3|. However the CI method supersedes the DS 
method in accuracy and so we will not consider the latter further here. We instead use scaling arguments to develop 
analytic solutions in the next section. 

IV. ANALYTIC SOLUTIONS 

In this section we employ analytic techniques to obtain solutions for the metric perturbations modes that are either 
far inside or outside the horizon on the branc. We begin in 3IVAI by deriving the behavior of perturbations in the 
quasistatic regime k/ Ha 3> 1, then we examine the opposite large scale limit k/ Ha <C 1 for cases where the expansion 
rate on the brane evolves as a simple power of the scale factor in flIVB| and finally in ^IVCfiVDl we turn to the 
behavior of superhorizon modes during dc Sitter expansion. 

A. Quasistatic Modes 

In this subsection, we examine the behavior of modes when kj Ha 3> 1 in the context of the quasistatic approximation 
and derive the first order correction to g in Ha/k. (First order corrections to the quasistatic limit of DGP have also 
been considered by Amin et al. |20|.) 

If we make the assumption in the master wave equation ||16a[ l that time derivatives of f2 can be neglected, we find 
that [ll| 

n w ci(l + zHy) k ' Ha + c 2 (l + eHy)- k / Ha (36) 

in the limit of kj Ha 3> 1. To ensure regularity of the master variable on the Cauchy horizon we need to set c\ = 
for the self-accelerating branch and c-i = for the normal branch. Therefore, the brane gradient becomes 



Ha 



'£ 



(37) 



This analysis justifies the quasistatic assumption in £|II El that the impact of R on the brane is negligible compared 
with {k/ Ha) 2 terms in the brane boundary equation (|17[) . 
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FIG. 2: Comparison of simulation results (black solid) to the quasistatic approximation at zeroth order (red long-dashed) and 
with first order Ha/k corrections (blue short-dashed). In both the self-accelerating and normal branches, the residuals Sg 
between the approximation and simulation results are roughly an order of magnitude smaller when first order corrections are 
included for a < 1. The first order correction becomes invalid as k/Ha — > as the mode exits the horizon during the de Sitter 
epoch o>l. Also note that we recover the general relativity result g — in the early time limit Hr c — > oo or, equivalently, 

Putting this into the expressions for the metric potentials (|18[) and expanding in inverse powers of k/Ha we obtain: 



$ + * 



K 2 pa 2 A 

k 2 
K 2 pa 2 A 

k 2 



nl 



7 4 (4iir c 74 - 373) Ha 



9Hr c j 3 
_ 74 Ha 
373 k 



O 



37l 
H 2 a 2 



k 



O 



H 2 a 2 



which in turn yields the first order correction to g from the brane gradient 



9 = 9QS 



e 74 [ 7l 2 + 3effr c (4ffr c7 4 - 3 73 )] He 
9ifr c 7| " k 



O 



H 2 a< 



(38a) 
(38b) 

(39) 



In Fig. [21 we show a comparison between simulation results and this formula when only the leading order term is 
retained (i.e. the standard quasistatic approximation) and when the next to leading order term is retained (i.e. a 
"first-order" quasistatic approximation) . 

The first order result gives a more accurate approximation to the simulations results in the k/Ha S> 1 regime. It 
extends either the accuracy or the regime of validity in k/Ha of the quasistatic approximation by about an order of 
magnitude. In the example in Fig.[2J it corrects the quasistatic approximation at the 10 -3 level for a < 1. 

For a>l, even modes that are far within the current horizon eventually exit the horizon. Note that the first order 
correction is unbounded as k/Ha — > 0. The numerical solutions on the self accelerating branch are also unbounded 
but the first order correction does not capture this behavior accurately. On the normal branch, the numerical solutions 
remain finite and again the behavior is not captured by the approximation. This problem applies equally to near 
horizon modes at a ~ 1. We therefore now turn to study the large scale behavior of g. 

B. Power- law Expansion 

In the supcrhorizon regime k/Ha <C 1, there is no universally-valid closed-form solution for g. In this subsection, 
we thus restrict ourselves to the situation where the expansion of the brane is well described by a simple power law 
in the scale factor a. We can then take 



H 2 ex a 



-3(1+10) 



H 



-°-(l + W )H\ 



(40) 



Here, the constant w is an effective equation of state parameter not to be confused with the equation of state of the 
matter on the brane. The scaling ansatz ([55)1 reduces the master wave equation ([TB]) down to an ODE for the scaling 
function G of the form [see H Eq. (43)] 

+ *>(*) +C0OG-O, (41) 



where x = eHy. The past Cauchy horizon in the bulk H- corresponds to 

2-hor — €-£1 2/hor 



l/(2 + 3w;), e = +l, u> > -2/3 

undefined, e = +1, w < -2/3 (42) 

-1, e = -l, 



while the brane is always at x = 0. The fact that the horizon position is undefined in the self-accelerating case when 
if < —2/3 is a consequence of the behavior of the Gaussian normal coordinates in the bulk. This is illustrated in 
Fig. [3l where we see that if the effective equation of state w < —2/3 at early times, t = constant hypcrsurfaces fail 
to intersect H_ when e = +1. (This also happens in the normal branch when w = —1, but this pathology manifests 
itself differently below.) We now seek solutions to Eq. (|41|l in the superhorizon k/Ha — > limit for x <E (0, x\ lol -) on 
the self-accelerating branch (when Xh or is finite) and for x £ (xhor, 0) = (— 1, 0) on the normal branch. 

A straightforward analysis of the poles of V and Q reveals that x = 1/(2 + 3w) and x = 2/(3w + 1) arc regular 
singular points of the ODE (|4"T| . The first of these is coincident with the bulk horizon when e = +1 and w > —2/3. 
The second singularity occurs outside the range x G [—1, 1/(2 + 3u>)] for all w, and hence is not relevant to our 
problem. 

In addition to the regular singular points, there is an irregular singular point at x = — 1; i.e., coincident with the 
bulk horizon when e = — 1. The nature of this singularity is apparent if we expand Q in a Laurent series about 
x= -1: 

^r \ 3 fc 2 1 + w ,.„. 

Q(x) = -- +••• 43 

AH 2 a 2 (1 + x) 3 

We see that this third order pole disappears if either w = — 1 or k/Ha = 0. 

We use the technique of matched asymptotic expansion to obtain the solution in the k/Ha -C 1 limit. Away from 
the pole at x = — 1, we can set k = in Eq. (|4T|) to obtain the outer (or near-brane) solution valid for w ^ —2/3 

G{x) = ci(l + xf + c 2 (l + x) 3/2 [l - (2 + 3w)x} {2p - 3)/{4+6w K (44) 

Here, c\ and ci are arbitrary constants. 

The outer solution is sufficient for the consideration of the self accelerating branch where x £ (0, Xhor)- When 
w > —2/3, the condition that £l(H-) = is equivalent to setting G(l/(2 + 3w)) = 0, which implies that c\ — when 
p > 3/2. Under these conditions, the brane gradient is 

For p < 3/2, strictly speaking no solution exists for constant w. However in the real universe, one would expect 
deviations from other components that break the w =const. assumption to make the c 2 solution finite at the horizon 
and c\ dominant at the brane. Hence we can assign R = p to the p < 3/2 modes. 

For w < —2/3 as in the case of the de Sitter regime considered below, a^or in principle diverges. In practice Xhor 
always remains finite at finite time given the preceding epochs of radiation and matter domination but increases 
without bound. 

Turning our attention to the normal branch, we require an inner (or near-horizon) solution close to the irregular 
singular point at x = —1. The Laurent expansion (|43[) implies that we cannot set k — a priori] i.e., the k/Ha 
contributions to the ODE arc not subdominant near x — — 1. To overcome this, we transform to the variable 



Ha / 1 + x ..„. 

^-k^w^r (46) 

In terms of this new coordinate, the horizon is at £ = and the brane is at £ 3> 1. The inner solution is specified by 
setting fc = in the transformed equation for £ and takes the form 



g(x(o) = e +3/2 



C3h/2-p ( 7 ) + c 4 K p-3/2 f 7 



(47) 
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(d) e = — 1, nonzero brane tension and 
no brane matter (de Sitter brane) 



(c) e = — 1, nonzero brane tension and 
brane matter with P = (CDM) 



(f) e = — 1, zero brane tension and 
brane matter with P = — |p 



FIG. 3: Conformal diagrams showing the behavior of the t — constant hypersurfaces in the space spanned by the (u, v) 
coordinates (p and P are the density and pressure of brane matter, respectively). The scaling ansatz is expected to be 
informative when these surfaces intersect the brane's past Cauchy horizon H-. This is indeed the case in the self-accelerating 
and normal branches when the brane transitions from early matter domination to a late time de Sitter phase, as in panels (b) 
and (e) respectively. However, when the brane undergoes purely de Sitter like expansion, as in (a) and (d), the t = constant 
curves fail to pierce T-L-. In these scenarios, it is impossible to impose the initial condition fi(%_) = in the context of the 
scaling ansatz. Finally, there are some choices of brane matter and tension where the behavior of the t = constant surfaces are 



qualitatively different on the two branches. An example of this is in panels (c) and (f), where the brane matter has P = 
In the normal branch t — constant surfaces cover H. 
is true for the self-accelerating case. 



and we expect the scaling calculation to be applicable, while the opposite 



where I and K are modified Besscl functions of the first and second kind, respectively. 

On the normal branch, the boundary condition that £l(7-LJ) = implies G = at £ = 0. Taking note of the 
asymptotic expansions of the Besscl functions for £ -c 1, 



h/2-p 



M i 



K 



P-3/2 



cxp 



(48) 



we see that we must set c^ — in the near- horizon G solution (|47[) . 

Finally we match the c\ inner solution (|4T|) to the outer solutions (|4"4"]) at £ 3> 1 but (1 + x) <C 1 to obtain the 
global approximate solution 



G( X (0) oc J^ +3/2 ^-3/ 2 (m 



(49) 



P > 3/2, 
\e +3/2 K p _ 3/2 (l/0[l - (2 + 3w)x}^P- 3 y( 4+6w \ V < 3/2. 

Substituting this global solution into the exact equation (I4T1) for G, we see that the residuals arc 0(k/Ha) and 
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converge uniformly in the interval < 1 + x < 1; moreover numerical integration of Eq. (|4Tj) confirms these solutions. 
Given £ 2 ex (1 + x), this solution implies that p > 3/2 matches onto the c x mode p < 3/2 the c 2 mode of Eq. (|4"4"|) 
for the normal branch. Before moving on, we note that these expressions for G in the normal branch only hold when 
the irregular singular point at x = —1 is present; i.e. when w > —1 and k is finite. In other words, the scaling 
approximation will not be informative in the purely de Sitter case w = —1, as could have been surmised from panel 
(d) in Fig. H 

To summarize, we have shown that under the scaling ansatz Q(t 7 y) = a p G(eHy) combined with the initial data 
G(eHyh OI ) = the brane gradient satisfies 



R = 



3- 


-P, 


e = +1, w > 


-2/3, p> 3/2, 


P, 




e = +1, w > 


-2/3, p< 3/2, 


P, 




f=-l,!ll/ 


-l,p>3/2, 


3- 


-P, 


e = — 1 , w 7^ - 


-l,p< 3/2 



(50) 



in the k/Ha — > limit. 

The preferred values for p are determined by the brane boundary equation through the density source. For example 
in the matter dominated phase when w — 0, the matter density fluctuation Aota acts as an external source to fib in 
the boundary equation (TIT)) so that for k/Ha <§C 1, the particular mode is p = 4. For the self accelerating branch the 
matter dominated solution is therefore R = — 1 and 



whereas on the normal branch R = A and 



3SH 



.9SH 



8Hr c - 1 



1 



2Hr c + 1 ' 
These relations close the perturbation equations on the brane. 



+1, 



e = -1. 



(51) 



(52) 



C. De Sitter Scaling 



The late time de Sitter phase carries an effective equation of state of W = — 1 and is a special case for the scaling 
solutions both on the self-accelerating and normal branch. The horizon on the self-accelerating branch goes to infinity 
while the horizon on the normal branch is no longer an irregular singular point of the master equation (|16a|) . 

Given the matter in the universe, a pure de Sitter expansion is never fully reached and it is important to consider 
the manner in which the de Sitter phase is approached. The Hubble parameter approaches a constant 



H* 
H 



y^A + ^r 



for a scale factor a 3> a* where H(a+) = \/2H- kr and 

Q m a- 3 = n A + (4 - 2V2)Q rc (l + ey/l + n A /Q rc ). 
(f2 m /f2 A ) 1/3 as usual. If Cl A /Q, r< , < 1 



Note that if f2A/^r c ^> 1, a* 



a* = 



a 



4(2 - >/2)n, 

i l rn 

(V2-i)n A 



1/3 



1/3 



-1, 



= -1. 



(53) 



(54) 



(55) 



In both cases Hr c also approaches a constant 



H*r c 




(56) 



We shall see that the residual effects from the preceding matter dominated phase determine the behavior of the 
perturbations on the brane. Likewise in Fig. [31 we see that the pathologies in panels (a) and (d) disappear with the 
addition of matter in (b) and (e). 



12 
Bulk solutions 



If we take a pure de Sitter limit where H = 0, the master equation (|16a[) becomes 



where recall ' = d/dlna. 

The general solution to the pure de Sitter master equation for k/Ha <C 1 is 

O = a p [ Cl (l + xf + c 2 (l + xf-P]. (58) 

This solution in fact corresponds to the w = —1 limit of the general constant w scaling solution (|44[) but in this case 
separability in a and x is guaranteed, not an ansatz. 

One might therefore expect the scaling results to hold: c\ = (R = 3—p) for the self-accelerating branch and c 2 = 
(R = p) for the normal branch for the fastest growing modes. However, direct application of the condition G(iEhor) = 
fails to be informative in both cases due to a change in the nature of the Xhor point. For the self-accelerating branch, 
its value diverges and for the normal branch it is no longer an irregular singular point of the master equation. The 
can also be seen in Fig. O where we see that if the brane undergoes a purely de Sitter expansion, the t = constant 
hypcrsurfaces do not intersect the past horizon and it is impossible to impose the initial condition O('H-) = in the 
scaling approximation. 

To restore the information lost in the pure de Sitter limit, we must consider the preceding matter dominated 
expansion as in panels (b) and (e) of Fig. [3] For the self-accelerating branch, the matter dominated phase leads to 
a finite but growing Xhor ~ aH^/Cl m H§. The c\ solution however does not depend on the evolution of H{a) and so 
even given this preceding phase its contribution to ft at the horizon increases as (1 + Xhor) p for p > 0. Hence the 
c\ contribution must be negligible near x ~ in order to satisfy the initial data G(xhor) = 0. The form of the c 2 
solution does depend on the evolution of H{a) and hence it is plausible that the preceding matter dominated epoch 
induces a correction of the form (1 — x/ x^ OT ) as it does in the matter dominated limit [see Eq. (|44[) ]. This correction 
has negligible impact near x = and so we again obtain 

i? = 3-p, e = +l, (59) 

for the solution on the self-accelerating branch. We shall see that this line of reasoning is borne out by the more 
detailed Green's function calculation below. 

On the normal branch, restoring a trace amount of matter has more dramatic effects. Defining 

h m = TTh—^- K a ~ 3 ' ( 60 ) 

Hi 

we obtain H' /H ps — (3/2)h m . In Eq. ©, the zero of n is shifted beyond the horizon x < — 1 and the master equation 
(|16ap regains an irregular singular point at finite k/Ha. 
We can transform to the variable 



Ha 1 + x 

e= -V^: (61) 

and again repeat the asymptotic matching of the interior k — and exterior k = solutions. As with constant w the 
interior solutions are given by Eq. (|47|) . For p > 3/2 they again match onto the c\ solution with R — p. For p < 3/2, 
it is important to consider an intermediate solution at finite h m but k = 0. The interior solution matches onto the 
(1 + x) 3 ' 2 rather than the (1 + x) piece of this solution and in general can stimulate both the C\ and the c 2 pieces of 
the exterior solution. From explicit solutions of these equations, we find that as h m — > the c 2 piece where R = 3 — p 
generally dominates given its larger growth with (1 + x). We shall see however that this leading order behavior is 
forbidden by the brane boundary equation. 

2. Brane boundary 

As in the matter dominated limit, certain values of the scaling index p are picked out by the boundary equation 
(f3"2"| . Once the matter contribution to the expansion becomes negligible n\p/ H 2 <C 1 and the matter conservation law 
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(|20|) requires that density perturbations obey A = Ao + Ai (a/a*) -2 where Ao and Ai are constants, i.e. to leading 
order, density fluctuations freeze out but with some stimulation of a decaying mode at the de Sitter transition. 

While this behavior of A is the same as in General Relativity with dark energy, the particular mode here does not 
drive the leading order behavior of the metric perturbations $ and iff. Note that the boundary equation ([T7]) and 
master equation (|16a[) requires a constant response in J7b 

2eHr cK lpa 3 
° b p^2 A o> ( 62 ) 

where the specific coefficient applies to R = p = 0. By virtue of Eq. (|T5)) . this mode caries no source to $ and *$> 
and as we shall see is the relevant particular mode for the normal branch [see Eq. (|83[) ]. Hence their evolution is 
typically governed instead by the homogeneous A source free modes in ilb- If the source free modes decay faster than 
p = —2, then the next to leading order scaling of Ax will dominate instead. On the self accelerating branch, we shall 
see that the fastest growing modes have p > and again the leading order particular contribution does not dominate 
the metric evolution. 

Now let us consider these homogeneous modes. We obtain from Eq. 



p(p 3) + -L_ + 2eHrc ^ = 0. (63) 

en 7'c en r c 



For the self-accelerating branch, the condition R = 3 — p (ci =0) gives 

p = 2,3--!- 1 e = +l, (64) 

Hr c 

where note that if fl\ = 0, Hr c —> 1 and both solutions return p — 2. This analysis justifies the selection of the p = 2 
mode in |3j. 

For the normal branch, taking R = p (c2 = 0) gives 

P = l,—^-, e = -l (65) 

n T c 

Here p < 3/2 and so in principle the R = 3 — p mode should dominate leading to a contradiction. If we instead 
assume R = 3 — p, Eq. (f6"3")l requires that either p = 2 or p = 3+ 1/Hr c ; in either case, p > 3/2 and we again reach a 
contradiction. Therefore, no solution of this type can satisfy both G(cchor) = and the brane boundary condition. 

Since the details of the matter to de Sitter transition are important for the R = 3 — p mode whereas they are not 
for the R — p mode, we conjecture that solutions with p < 3/2 should exist if corrections near the singularity prevent 
the C2 mode in Eq. (|58[l from completely dominating. 

If we further assume that there exist modes where R = p is the dominant mode at p < 3/2 we obtain the modes 
in Eq. (|65j). In principle though, we can obtain a spectrum of modes with mixed R = p and R = 3 — p behavior. 
We shall see that the numerical CI solution implies that p = 1 is indeed the dominant mode for Q. However, this 
mode has no effect on the brane metric parameters $ and \P. Hence the important k/ Ha — > mode is p = —l/Hr c , 
but that is a decaying mode in $7. This opens up the possibility that fc-dependent modes eventually dominate the 
evolution in the de Sitter epoch. To study this behavior and verify the conjectures involving the c-i mode for both the 
self- accelerating and normal branch, we require a more complete Green's function analysis. 

D. De Sitter Green's Function 

To gain further insight into solutions in the de Sitter epoch, we can recast the perturbation equations as a canonical 
scattering problem via the transformation 

n(t, y) = (1 + eHy) 3 / 2 a 3 / 2 (tMt, y), (66) 



and define a new bulk variable z: 



ln(l + egy) e egz 

Z= eH ' y= ^lF- (67) 
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The bulk manifold is defined by z £ [0, oo) for both branches. The new bulk variable satisfies a simple wave equation 
and boundary condition 



d 2 d 2 k 2s 



dip 

dz 



cr. 



9V 
2Hr c - e V dz 2 



eH{3Hr c - 2e) 
4(2Hr c - e) 



3/2 



-<£>b 



2eK4r c pa 

k 2 {2Hr c -e) 



A. 



(68a) 
(68b) 



where we have set <fb{t) — <p(t, 0) and made use of Eq. (|68al) to remove time derivatives of tp. The equation of motion 
for A reduces to 



2HA 



n\pHr c 
2Hr r - 



A 



efc 4 



3(2Hr c - e)a 7 / 2 



¥>b- 



(69) 



Examining equations (|68[) . we see that A can be viewed as a brane- localized source for the bulk ip field. From linearity, 
it follows that the general solution for ip should be of the form 



ip = <p {h) + ip {p \ 



(70) 



where tp^ is the general homogeneous solution for Eq. (|68[) when A = 0, and ip^ is a particular solution when A is 
nonzero and determined from Eq. (|69[) . 

Let us now concentrate on the solution of the homogeneous system where ip^ h ' and its first time derivative ip( h > are 
known at an initial time, say t = 0. Then, the value of ip( h > at some later time t and position z is given by 



<pW{t,z) 



dz' 



i P ^(i/,z')-^-G_(t,t'-,z,z')-G-(t,t';z,z')^-<p^(i/,z') 



dt 
er, 



or 



t'=0 



2Hr c 



^\t',z')^G-{t,t l ^z,z')-G-{t,t , -,z,z')^ { - h \t',z') 



dt' 



dt 1 



(71) 



t'=0,z'=0 



where G_ is the retarded Green's function. Though this formal solution of the initial value problem only involves G_, 
we will also consider the advanced Green's function G+ to facilitate comparison with the scaling solutions of ijIVCI 
Both Green's functions satisfy 



W 2 ~& 2 + V 2 ^ 



with boundary conditions 
= 



d er c d 2 eH(3Hr c - 2e) 

~dz~l + 2Hr c - e ~dz\ + 4(2iJr c - e) 



G±(t,t';z,z') 



z<=0 



< Az ± At 



G ± (t,t';z,z') = 0, 



(72) 

(73a) 
(73b) 



where z < = min(z,z'), z> = max(z,z'), At = t — t' , and Az = \z — z'\. The boundary condition (|73bl) ensures that 
the initial condition S1("H_) = will be satisfied by Eq. (|7Tj) for t > 0. Note that T-L- corresponds to t — > — oo and 
z — > oo, yet the boundary condition (|73b[) is imposed at finite values of the coordinates. This is in contrast to the 
scaling approach, where the condition Q = is enforced on T-L- directly. Therefore, unlike the scaling solution, the 
Green's function approach is not compromised by the fact that Gaussian normal coordinates do not cover the past 
horizon for a de Sitter brane. That is, it is not necessary to impose conditions on the Green's function at 7i- to 
ensure that the solution (f7T|) is consistent with Vt{T-L-) =0 . 

Using separation of variables and other standard techniques, the advanced and retarded Green's functions can be 
expressed as 



G±(t,t ; z, z) 



+oo 



dvG±(t, t';z, z'\v), 



g±(t,t';z,z'\v) 



T v (t)TW)[f T (v) 



2ivHfr(y) 



f ± {v)e ±lvH ^ z+z ">] 



(74a) 
(74b) 
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Green's function 


resonant mode 


frequency v 


temporal scaling p 


brane gradient i? 


retarded G_ 


A 2 
B m {m=l,2...) 


-it/2 
-i(3e/2 - u) 

im 


(3 + e)/2 

3(l + e)/2-w 

3/2 -m 


1 

7i(9-4eoj-4m 2 )/8 


advanced G + 


Ci 

c 2 

An (to = 1,2...) 


it/2 
i(3e/2 - o>) 


(3 - e)/2 
3(l-e)/2+w 
3/2 -m 


1 

7i(9-4eoj-4m 2 )/8 



TABLE I: Resonances of the retarded and advanced Green's functions as determined by the poles of Q- and Q+, respectively. 
(to = 1/Hr c ). Note that the A\^ and Ci,2 excitations give rise to late time resonant modes of the form fi oc a p (l + eHy) R , 
while the B m and D m modes have more complicated y dependence. 



The two terms represent propagation directly from z' to z or through reflection off of the boundary. The functions 
T v satisfy 



S,+K+h^ 



dvT v (t)T:{t') = 5{t-t l ), 



and f± follow from the boundary condition (|73a|) : 

f±(u) = (2ev ± i)(3iHr c - 2ie ± 2eHr c v). 
Explicitly, the temporal mode functions arc given by 

1/2 



k V" / Hv V'" ( k 

2~H) \&mh.'Kv) ~ iv \a~H 



(75) 



(76) 



(77) 



Here, ./_«, is the Bessel function of the first kind with order —iv. 

Now, in order to evaluate the righthand side of Eq. (|74a[) . we close the contour of integration in the upper- half of 
the complex v plane when =pAi > Az, and in the lower-half plane otherwise. This gives 



G±(t,t'-z,z') = e(TAt-Az) 



2™ V Res £± (£,£'; z,z» - / dvG±(t,t'; 



z,z'\v) 



(78) 



In this expression, the sum is over the residues of the poles of G±(t,t';z,z'\v) in the upper-half of the complex v 
plane, and the contour T + represents the infinitely large semi-circle used to close the integration path. The poles of 
G±{t, t'\ z, z'\v) represent resonant excitations that will dominate the late time behavior of our system, and are listed 
in Table HI The A and C resonances listed in the table arise from zeros in f±{v) in Eq. (|76|) . while the B and D 
modes are associated with divergences in the T v (t) functions. In what follows, we will neglect the T + integration and 
represent the Green's function as a sum over resonances. 

It is instructive to relate the A and C resonances to the preferred modes in the scaling approach. The scaling 
solution of Eq. ([58| corresponds to 



<p(t,y) - Cie (P-3/2)H(t+^) +C2e (p-3/2)H(t-ez)^ 



(79) 



Each term represents a wave traveling towards or away from the brane depending on the choice of branch, and the 
boundary equation satisfied by ip can be interpreted as a reflection condition on these waves. Comparing the resonances 
from the poles of G± to Eq. (fT9")l . we see that the A 12 resonances of the retarded Green's function correspond to scaling 
solutions with c\ = in the self-accelerating branch and ci = in the normal branch. Conversely, the C\^ advanced 
resonances have C2 = and c\ = for the self-accelerating and normal branches, respectively. The causality condition 
(|21j) picks out the retarded solution in each branch, which in turn implies that c% = and R = 3 — p for the self- 
accelerating branch, while C2 = and R = p for the normal branch. This analysis verifies our conjectures of the 
previous section based on extending the dc Sitter scaling solutions to include trace amounts of matter to specify the 
initial data. 

In order to facilitate the comparison of simulation results with the predictions of the Green's function analysis, we 
hereafter restrict our attention to the retarded Green's function and assume that the field point (t, z) is on the brane 
and deep in the de Sitter era; i.e. k/Ha(i) <C 1. The asymptotic form of G_ in this limit is 



G-(t,t';0,z') 



8ieHr c 

7l7T 



+ OG 



dv 



2Ha{t) 



^ r O) -iuHz' T 

/+0) 



Ha(t') 



Ha{t) 



<1. 



(80) 



16 



This is the form of the Green's function that one would use to determine the behavior of £7b at late times. It is also 
interesting to consider the dynamics of modes which are in the superhorizon regime at the beginning of the de Sitter 
era: k/Ha(t') <C 1. These are the modes which never actually enter the Hubble horizon. It is possible to explicitly 
evaluate the residues of Q- in this case to obtain 



G~{t, t'; 0, z') » ©(At - Az) I A x 



a (t) -Hz 

o(f) 



72 



Ao 



a>(t) & -Hz 



a(f) 



3e/2-l/Hr c 



£*« 



k 2 e Hz ' 
H 2 a(t)a(t') 



(81) 



where A\, A2 and B m are functions of Hr c . One can easily derive a similar expression for the advanced Green's 
function with coefficients C\ , Ci and D m corresponding to the advanced modes in Table HI 

The Green's function ([8lj) reveals how scale-dependence outside the horizon can arise. The A\ and Ai contributions 
to G- arc independent of scale k/Ha. This is consistent with what we would expect in General Relativity for a mode 
that has k/Ha < 1 at the beginning of a de Sitter era: the evolution of the mode (as mediated by the Green's 
function) is expected to be independent of k. On the other hand, the contributions from the B m modes carry an 
explicit /c-dcpcndcncc (this arises from the fact that Ji V and J_i„ are linearly dependent when iv is an integer). 
This (k/Ha) 2 suppression is also familiar from General Relativity and represents the leading order correction to the 
gradient approximation. However, on the normal branch, when combined with more slowly growing A\ and Ai modes, 
the net growth rate of the metric fluctuations for superhorizon perturbations becomes fc-dependent unlike in General 
Relativity. 

We now turn our attention to finding a particular solution to the inhomogeneous equation (|68[) with the A source. 
In principle, one could use the retarded Green's function G_ to find ip^- p \ but it is easier to just solve for it directly 
in Eq. (|68[) by making a physically- motivated ansatz. In this calculation, we work in the a — > 00 limit and retain only 
leading order terms. We make the following "outgoing-wave" ansatz for the bulk field: 



<pW oze wH{t - z \ 
Substituting this into Eq. (|68|) and (|69|) . we find that there are two values of v consistent with A^O: 

v = 3z/2 or 7i/2. 



(82) 
(83) 



The two distinct values of v give rise to solutions for A and D, h 
A(i)=A + Ai(a/<T 2 



(p) 



Ar 



la 



A a 



(p) _ „3/2,>) _ r c n A pa 



2 pa 3 
Hk 2 



3ifr c - 

-2A 



1 2a 2 (5Hr c -l)' 
2al Ai 

'3a 2 (2Hr c -1)' 



e = +l, 
e = -l, 



(84) 



which confirms the scaling expectation [see Eq. (|62|) ]. 

Having now determined how to obtain both homogeneous (p^ and particular solutions (p^ of Eq. (|6"5|) at late 
times, we can write down the following formulae for various quantities of cosmological interest 



fit 



I>i" f 







(p) 

b ' 



$«y$ (!) + $ (p) , jsiV^'i+f 



( P ) 



(85) 



Here, Pi is the temporal scaling index of the i th resonance of the retarded Green's function (i.e. the p- values for the A\ : 2 

and B m modes listed in Tabled, ^0 arc constants with dimension (length) 2 related to the amplitudes (A\,A2,B m ) 
and recall a* is the reference epoch at the beginning of the de Sitter phase [see Eq. (|54[)]. The resonant contributions 
to the metric potentials and comoving curvature perturbation are 



$(i) = niVfe - 1) 



2a*(2#r c e-l) 



*W = ( p . t _ i)$W 



For Pl £ 1, 



l + O 



(») 

ffsH 



Pi-1 



H 2 a 2 



Pi 



l + O 



/ k 2 



\H 2 c 



(86a) 
(86b) 

(87) 
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branch 


matter era 


Ax 


^2 


Bj 


particular mode 


e=+l 
e= -1 


+9/(8Hr c - 1) 
-l/(2Hr c + 1) 


oo 
n/a 


3^r c - 1 
Hr c -1 
-l/(2Hr c + l) 


1/3 
1/3 




-l/(2Hr c + 1) 



TABLE II: The value of the quantity g in the superhorizon regime k/ Ha <C 1 in either branch. 

for each mode on either branch. This scaling is in fact a direct consequence of £' = in Eq. (|28|) for pi ^ —2. 

On the self-accelerating branch with tt\ = 0, the growth index p t = 2 for both A\ and A 2 . This large growth 
rate makes them dominates over B and particular modes. In fact <?sh ~^ °°j which explains its divergent behavior in 
Fig. [21 For the normal branch, the Ai mode has pi = 1 and thus no contribution to either $ or "J. The fc-independent 
solution that dominates is then A 2 but its growth index pi = —1/Hr c is smaller than the B\ mode pi = 1/2. Thus 
the fc-dependent B\ mode dominates at late times for any finite k. 

On the normal branch even the particular mode can be important. The particular contribution to the various 
quantities are given by 

3 An 



--^^l\Jl^ 



k 2 "2/ 



k 



e - 


= +1, 


a*\ 2 2Hr c A 1 
a) AH 2 r 2 c -V 


= -1. 




r = +l, 


a+\ 2 2(Hr c + l)A 1 
a ) AH 2 r 2 - 1 ' 


e= -1. 



3 \HaJ 2Hr r 
3 Ap 
(p) _ Kipa'Hr c I '2 3Fr c -l' , . 




3 V^ a / 2Hr c + l 

Note that for the normal branch the leading order An term vanishes and contributions are suppressed by (k/Ha) 2 . 
Hence for modes that are outside of the horizon at the de Sitter transition (k/Ha+) 2 < 1 we expect that the main 
contribution from the particular mode comes from Ai. Under these assumptions the leading order scalings give 

e = +l 

< 89 > 

For the normal branch, the growth index p( p > = —2 and hence can become the leading order /c-independent term 
if —1/Hr c < —2. However the value of <?sh remains the same as for the A 2 mode and is in fact the same as the 
matter dominated scaling. Therefore on the normal branch, the leading order fc-independent term can be modeled 
a s ,9sh = —l/(2Hr c + 1) regardless of which of those modes actually dominates the solution. This simplification was 
employed by [9( to model perturbation evolution to the present epoch and obtain observational constraints on Hor c . 
In Table HU we list the values of c/sh associated with each of the modes appearing in Eqs. (|55]l and (JSSJ), as well as 
the values of <7sh in matter domination from JjIVBI These relations supplemented by the quasistatic analysis complete 
the analytic description. 

V. FITTING FUNCTIONS 

Given the analytic description of perturbations in the matter and de Sitter epochs above and below the horizon in 
QYVi we now devise global fitting functions for the perturbation evolution that bridge the transition and apply across 
all scales and epochs. We focus primarily on the normal branch but also test existing fits for the self-accelerating 
branch that have been used to show its predictions are in conflict with CMB data p, [l5| . 

A. Normal Branch 

As discussed in £|II Ci once g = ($ + v I / )/( < ! > — *) is determined, solving for the metric itself is a simple matter of 
applying conservation of the comoving curvature on large scales or the modified Poisson equation on small scales. 
Moreover g = for General Relativity without anisotropic stress sources and an accurate quantification of its value 
in DGP is important for observational test of gravity. 
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FIG. 4: Metric ratio g for purely superhorizon modes for normal branch (e = —1). Analytic fit from Eq. f|93f) (dashed lines) 
match the numerical CI calculation well for a wide range of parameters. Note that even on superhorizon scales g is scale 
dependent due to the relative importance of the A, B, and particular modes. The current epoch is denoted by the arrow. 



We begin by describing the superhorizon behavior. As shown in the previous section, there are several modes of 
importance during and after the transition to a de Sitter expansion. Although the growth rate of all modes are 
independent of scale above the horizon, their initial amplitudes are not. In particular, after the de Sitter transition 
at a+ where H(a — > oo) = H+ = H(a+)/y/2 [sec Eqs. (|53|) and (|54)) ]. there are scale-free modes and modes that come 
from higher order terms in the gradient approximation that are suppressed by powers of (k/ H+a*). The leading order 
mode By has a faster growth rate and though initially suppressed will eventually dominate on all scales [see Eq. ([51])]. 
B\ requires g = 1/3 whereas for the scale free A resonant and density-driven particular modes g = —l/(2Hr c + 1). 

Given these relations and the growth rates of the three modes, we seek to describe the superhorizon behavior as 



9sn(a, k) 



{k/H^f - [^(q/q,)- 1 / 2 - 1 /^)^ +K 2 (a/a^/ 2 ] 
3(fc/^q,) 2 + [2H{a)r c + l][jr 1 (o/a*)- 1 /a-V*(«)--c + K 2 (a/a*)-^ 



With 



K x = 13.2 



#*r c + 1 



ff*r c 
K 2 = 30.9(2fl> c + l) 2 , 



(90) 

(91a) 
(91b) 



we are able to fit the results of the CI numerical solution very well from matter domination through dc Sitter expansion 
for 240 models with 0.1 < H+r c < 10, 0.01 < fl\ < 1, and k/H±a± < 1. Examples of the fit are shown in Fig. [4j 
The subhorizon form of g approaches the quasistatic behavior at k/ Ha ^$> 1 



9Qs(a) 



31-2eHr c (l + H'/3H) 



(92) 



Note that in the dc Sitter epoch, modes exit the horizon. In choosing a functional form for mediating the transition 
it is useful to note that the first order predicted solution is g = <?qs(1 + 0(Ha/k)) from JjIV Al Unfortunately it does 
not suffice to simply match a first order correction to the superhorizon modes at horizon crossing. Instead we find 
that the following interpolation 



g{a,k) 



gSH + ffQsA" 3 

1 + K 3 



(93) 
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FIG. 5: Metric ratio g for modes that are currently horizon and sub-horizon scale for the normal branch (e = —1)- Analytic fit 
from Eq. (|93|l (dashed lines) match the qualitative features of the numerical CI calculation across the full range of parameters. 
In particular, the fit is designed to approach the correct behavior for kj Ha <g 1, k/Ha 3> 1 and a/a+ 3> 1. Note that currently 
subhorizon modes exit the horizon in the future de Sitter epoch. The current epoch is denoted by the arrow. 
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FIG. 6: Observable lensing potential ($ — ^0/2 on normal branch (e = — 1). Analytic fit to the metric ratio g compared with 
the numerical CI calculation (top panel). In spite of imperfections in the p-fit for currently horizon scale modes, fractional 
differences (bottom panel) are < 2% for a < 1. The current epoch is denoted by the arrow. 
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FIG. 7: Observable lensing potential ($ — \&)/2 on self-accelerating branch (e = +1). Analytic fit to the metric ratio g compared 
with the numerical CI calculation (top panel). Fractional differences in the observable potential are < 2% for a < 1 (bottom 
panel) . 



with 



I\\ = 0.18 



Ha 



<' r 



(8.1a/a,) E 



3 5 + (k/H+a+) 



^■(3H±r c + l) 



(94) 



better describes the k/ Ha ~ 1 horizon-crossing regime while providing a negligible mismatch in the k/ Ha 3> 1 limit. 
Note that for k/Ha 3> 1 and a/a+ ^> 1, the first order correction of Eq. (|39| gives (g — gQs)/gQS = — 2(i/+a/fc) 
whereas the fit gives (g — <7qs)/<7qs = — l-7(H*a/k)(l + H*r c )/(3 + H±r c ). That the functional form in k is correct 
and the H±r c corrections are bounded from 1/3 to 1 ensures that that we do not have a runaway mismatch for any 
set of parameters. 

As shown in Fig. [5j the fit captures the qualitative features of the numerical solution, including the asymptotic 
behavior at a/a* ^$> 1. Although imperfect, the fit suffices for observational constraint purposes. In Fig. El we show 
the performance of the fit for the metric combination involved in gravitational lensing and gravitational redshifts 
($ — 4 f )/2. Here we have set the transition scale between constant comoving curvature and Poisson-like behavior to 
c r = 0.15 [15(. Note that fractional errors arc below 2% for a < 1 for all scales. 



B. Self- Accelerating Branch 

Similar fitting functions for the self-accelerating branch were given in [15| for 
scales 



ffSH(a) 



9 



8Hr c 



1 



1 



0.51 



Hr c 



1. 



< 1 and Oa = 0. For superhorizon 



and the interpolation factor in Eq. (J93J) is given by K3 = (O.lAk/Ha) with the transition scale set by Cr = 1- In 
Fig. we compare the predicted lensing potential to the numerical integration up to the present epoch. Given that 
we have shown that gsn actually diverges deep in the de Sitter epoch (see Tab. |TT] and Fig. [5]), it is not surprising 
the errors in the fit increase toward the current epoch. Still those errors are < 2% at the time scales relevant for 
observational constraints. This should be compared with the much larger change in the potential of a factor of ~ 2. 
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VI. DISCUSSION 

We have provided analytic solutions for the linear evolution of metric perturbations in the DGP modified gravity 
models in various regimes where scaling and Green's function techniques are applicable. We have elucidated the 
nature of the coordinate singularities and initial data in the bulk as well as their effect on perturbation evolution on 
the brane. 

Interestingly, even on superhorizon scales, the evolution of metric perturbations is no longer necessarily scale free. 
In the late-time de Sitter phase on the normal branch, several resonant modes are excited with different growth rates 
and amplitudes. The epoch at which the fastest growing mode dominates depends on scale. On the other hand, 
for reasonable parameters where the de-Sitter phase has only recently been entered these complexities are mainly 
manifest in the future and do not affect observational tests [9| . 

Based on these analytic solutions, we have devised convenient fitting functions for the evolution that bridge the 
various spatial and temporal regimes for the normal branch. These forms are valid for essentially the whole range 
of parameter space and have been tested against numerical integration of 240 models spanning 0.1 < H+r c < 10, 
0.01 < 51a < 1 with multiple fe-modes for each. For the self-accelerating branch we have verified the accuracy of 
existing fitting formulae for the cases of interest [15[ . 

Compared with a direct numerical integration of the bulk equations, these forms are accurate at the percent level 
which is sufficient for current and upcoming observational tests of the DGP scenario from the cosmic microwave 
background and weak gravitational lensing. 
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